NASA Contractor Report 189636 
ICASE Report No. 92-14 






ICASE 


FAST WAVELET BASED ALGORITHMS FOR LINEAR 
EVOLUTION EQUATIONS 


Bjorn Engquist 
Stanley Osher 
Sifen Zhong 


Contract No. NAS 1-1 8605 
April 1992 

Institute for Computer Applications in Science and Engineering 
NASA Langley Research Center 
Hampton, Virginia 23665-5225 

Operated by the Universities Space Research Association 


NASA 

National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton, Virginia 23665-5225 

(NASA-CR-139636) FAST WAVELET 8ASE0 N92-24983 

ALGORITHMS FOR LINEAR EVOLUTION EQUATIONS 
Final Report (ICASG) 29 p CSCL 12A 


Unci as 

G3/64 0091243 



Fast Wavelet Based Algorithms for Linear 
Evolution Equations 

Bjorn Engquist, Stanley Osher 1 , and Sifen Zhong 
Department of Mathematics 
University of California, Los Angeles 
Los Angeles, CA 90024 


Abstract 

We devise a class of fast wavelet based algorithms for linear evolution equations whose 
coefficients are time independent. The method draws on the work of Beylkin, Coifman, and 
Rokhlin [1] which they applied to general Calderon- Zygmund type integral operators. We 
apply a modification of their idea to linear hyperbolic and parabolic equations, with spatially 
varying coefficients. A significant speedup over standard methods is obtained when applied 
to hyperbolic equations in one space dimension and parabolic equations in multidimensions. 
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1. Introduction 


During the last few years a number of fast computational algorithms have been developed 
for elliptic problems. These are techniques for which the number of arithmetic operations 
needed are close to linear as a function of the number of unknowns. Examples of algorithms 
of such complexity are multigrid methods and the so-called fast Poisson solvers. The fast 
multipole method and wavelet based methods for elliptic problems formulated as integral 
equations alsoTielong to this category [8], [1]. 

There has not been the same progress for hyperbolic and parabolic methods. In general 
classical numerical techniques for these problems are already optimal. 

Consider a system of evolution equations. 

d t u + L(x, d x )u = f(x), x 6 0 C R d , t > 0, 

( 1 . 1 ) 

u(x, 0) = u 0 (x), 

with boundary conditions, where L is a differential operator. 

An explicit discretization of this problem typically takes the form, 

u] « u(xj,t n ), t n = nAt, 

Xj — (ji Axj, . . . , jdAxd) 


(1.2) u n+1 = Au n + F, 

u° = u 0 , 

u,Fe R N “, At = const. | Ax| r . 

The vector u n contains all the unknowns u" at time level t n . For simplicity we shall assume 
j v = 1,2, . . . ,N in all dimensions v = 1, . . . , d. 

The matrix A is ( N d x N d ) with the number of elements ^ 0 in each row and each 
column bounded by a constant. Every time step requires 0{N d ) arithmetic operations and 
the overall complexity for a time interval of 0(1) is of the same order as the number of 
unknowns, 0(N d+T ). 

There are, however, some fast methods based on the analytic form of the solution opera- 
tor. In [3] the multidimensional heat operator, with u 0 and / both zero, but with inhomoge- 
neous boundary data given at M points, was treated. There the closed form of the solution 
evaluated at M points at time level N was obtained in O(NM) rather than 0(N 2 M 2 ) op- 
erations. Also, in [4], the same authors obtained an algorithm for evaluating the sum of 
N Gaussians at M arbitrarily distributed points in 0(N + M ) operations. So far, their 
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interesting method appears to need an explicit analytic representation of the heat kernel, 
effectively ruling out variable coefficient problems. 

The formula (1.2) has a simple closed form solution 

71 — 1 

(1.3) u n = A n u 0 + AVf - 

u=0 

This form can be used to compute the solution A n uo, for F = 0, in log n steps, (n = 2 m , m in- 
teger; here and throughout, log n = log 2 n) by repeated squaring of A : A, A 2 , A 4 , A 8 , . . . , A 2 "*. 

Unfortunately the later squarings involve almost dense matrices and the overall complex- 
ity is 0(N 3d log N) which is larger than that using (1.2) directly. 

For an appropriate representation of A in a wavelet basis all of the powers A" may be 
approximated by sparse matrices and the algorithm using repeated squaring should then be 
advantageous. 

We shall consider the following algorithms for the computation of the closed form solution 
(1.3) of the inhomogeneous problem in m = logn steps, 

B := SAS~ l 


(1.4) 


C I 

C :=TRUNC(C + BC,e) 
B := TRUNC(BB,e) 


(iterate m steps) 


u n := S~ 1 (BSu° + CSF). 

The matrix S corresponds to a fast transform of wavelet type and the truncation operator 
sets elements in a matrix to zero if their absolute value is below a given threshold. 

(1.5) A = TRUNC(A,e) : j * 

It is easy to see that algorithm (1.4) is equivalent to (1.3) for e = 0. This is not so for e > 0 
and also for F ^ 0. We shall however show that it is possible to choose £ small enough for 
the result of (1.4) to be arbitrarily close to (1.3) but still with very few arithmetic operations. 

For a fixed predetermined accuracy level the computational complexity to calculate a one 
dimensional hyperbolic equation can be reduced from the standard 0(N 2 ) to 0(./V(log N) 3 ). 
The extra cost per time step is minimal. This also makes it possible, as a curiosity, to use 
algorithms which are unstable in the traditional sense. 

Our technique is even more favorable for parabolic problems. A d-dimensional explicit 
calculation with standard complexity 0(N d+ 2 ) may be reduced to 0(N d (\ogN) 3 ). 
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The algorithm (1.4) can be extended to some problems with time dependent data. In this 
case, we clearly need to compress the information in the data such that not all the 0(N ) 

values in, e.g. the inhomogeneous term f(xj,t n ) are needed. 

One simple but important application of this type is from optics or electro-magnetic 
scattering with a time periodic source. If k points are needed to resolve one time period, we 
can group k time steps together 

(1.6a) u n+k = A k u n + J2 A j F n+k +j- 1 . 

3=0 

where 


(1.6b) 


/; = Af /(<«). 


This equation is now of the type (1.2) with time step kAt and with inhomogeneous term 


k - 1 

(1.6c) F = ^2 ^ Fn+k+j-l- 

3=0 

In sections 2 and 3 we shall discuss the analytical properties of the algorithm. Numerical 
examples are presented in section 4. 


2. Hyperbolic Problems 

Consider first the simple one dimensional scalar advection equation, 

d t u 4 -ad x u = 0, a > 0 

( 2 . 1 ) 

u(x,0) = u 0 (x), 0 < x < 1. 

The functions uq and thus u are assumed to be 1-periodic in x. The solution of (2.1) is given 


by: 


(2.2) u(x,t) = uo(x - at). 

The different rows of A v in a numerical solution of (2.10) will represent approximations of 
the Green’s function G below, 

/ OO 

G(x,y,t)u 0 (y)dy, 

-OO 


(2.3) 


u(x,t) 



X - y - at)u 0 (y)dy. 
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Let y>j be a truncated wavelet expansion of a (5-function with an orthonormal set of compactly 
supported wavelets, 

8{x) ~ <pj(x) = ^OL 3 kV]k{x) 

<pjk{x) = 2~^(2~ 3 x - k + 1) 

The choices of ip(x) will be discussed below. Assume that the rows of A v are discrete 8- 
functions, i.e. just one element is nonzero and large. For each level j = 1,2 there are 
only a finite number of ajk ^ 0. With J = m = log N there is only log N of all ajk ^ 0. 

Thus each row in B, (1.4), has log N elements, bjk ^ 0. The matrix B 2 is also a transform 

of an idealized matrix A v and will have N log N elements different from zero. This means 
that each iteration step in the algorithm (1.4) produces 0(iV(log N) 2 ) flops when F = 0. 
We have assumed that calculations are only carried out for those B 2 elements which are 
different from zero. In practice a slightly larger number of elements needs to be computed 
and then truncated. This corresponds to the case when the location of the <5-functions is 
only approximately known. Compare the wavelet technique for Burgers’ equation by Maday, 
Perrier, and Ravel [6]. 

Each row of C, (1.4), is a transform of a step function, 

\ f const. 0 < x < at, 

C <*> = { 0, else _ 

This function can also be represented by log N wavelets and thus the overall cost is 
0(N(log A0 3 ). 

In numerical computations the rows of A " are only approximations of (5-functions. If an 
upwind scheme, 

uf 1 = Uj - A(u" — Uy.j), 

(2.4) u° = u 0 (xj), j = 1,2, ... AT, 

A = aAt/Ax < 1, 

is used A will have the form, 

1 — A 0 
A 1 - A 0 

0 A 1 - A 0 

0 0 A 

The matrix A v will have Toeplitz structure. Each row is still an approximation of a <5- 
function. The first order smoothing effect of (2.4) is given by the modified equation, [5], 

(2.5) d t u + ad x u = (aAx/2)^u. 



A 

0 

0 


1 - A 
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Equation (2.5) is parabolic with a fundamental solution of the form, 


(2.6) G(x — y,t) = (2iraAxt) * exp(— (x — y — at) 2 j{2aAxt)). 

Compare the solution formula for parabolic problems (3.2). 

Each row of A u is thus a close approximation to the function G(x — y,t) above. The 
computational complexity of the algorithm (1.4) depends on how many wavelets are needed 
to represent G(x — y,t) as a function of x, (0 < t < T) with a given accuracy. 

Higher order accurate (say order 2p-l) dissipative finite difference approximations to (2.1) 
are usually modelled by the equation 


with k p > 8 > 0,(5, independent of Ax. 

The fundamental solution for this parabolic equation is: 

1 f°° 

G p (x , t) = — / d£ exp (i£(x - at) - fc p (Ax) 2p ~ 1 £ 2p f)). 

Z7T J — oo 

The key estimate we shall obtain here (and which we certainly do not claim is new) is: 

(± 

\dx 

uniformly in 0 < t and Ax and for all nonnegative integers m. 


(2.8) 


X 


G p (x + at,t) 


<c, 


^m,p 


(2.7) 


u e + au x = (-iy +1 kf,(Ax) 2p 1 (A 


Proof of 2.8. We wish to bound 


— r {iO m x m+x e^ x ~ k ^ x ^ hev d(. 

27T J -oo 

_ / n \ tti + 1 

= _L r ( 1) [ e 

2tt 7-00 U^j r 

= ir 

2tt J - 


m e -k p (&x) 2 r-H 2p t 


d( 


exp 


i£x 


XtiAxyp-'kpy/ip 


)] [dP^i 


d(. 


The result is now clear. Also, an inspection of the right hand side of the above shows that 
G m ,p can be chosen to be arbitrarily small if f(Ax) 2p_1 is large enough. 


Remark Rl. Let the general space dependent coefficient, one dimensional system of hy- 
perbolic equations 

u t + A(x)u* = C(x)u , 
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where u is an £ vector, A is a uniformly diagonalizable smooth £ x £ matrix, with all real 
eigenvalues A, (a:), and C(x ) is smooth, be approximated by a dissipative finite difference 
scheme of order 2 p — 1. Typically, its model equation is a systems version of (2.1) 

u x A(x)u x = C(x)u + (-l) p+1 (Ax) 2p_1 P ^x, j u 

where (— l) p+1 P(x, is a 2p order elliptic operator. A more involved argument shows that 
the fundamental solution satisfies an estimate of the type (2.8) with the expression x + at 
replaced appropriately by solutions of ^ = A ,(x) x(0) = x, i — !,...,£ and with C m , P 

possibly growing in time like C m<p e kt for k fixed. 

Our numerical procedure involves the compression of the matrix A u , which for the purpose 
of analysis only, we shall view as the discretization of the fundamental solution for either 
(2.5) or (2.7), 

( A n ) jk = G(xj,y k ,t n ) 
where the interval [0, 1] is discretized via 

j 


Xj 


N' 


j — 1, . . . , A, N = 2* 


[0, 1] x [0, 1] is discretized via (xj, y k ), and t n = nAt = nXAx , n = 0, 1, . . . . 

We now adapt the terminology, notation, and results of [1] to this unsteady problem 

( 1 . 1 ). 

Finite difference schemes approximating (1.1), e.g. (2.4) are regarded as acting on a 
vector which is to be viewed as approximating u(x,0) on the finest scale: 


s° = 2^ J <p(2 1/ x — k + l)u(x,0)da 
= J f(x)ip l/k {x)dx. 


All functions, both continuous and discrete, are extended periodically: 

u(x, t ) = u(x + 1, t) 


etc. 


The function satisfies 


<t° = q° 

S k+N — s k 


2m— 1 

<f{x) = Y, h p+ iip{2x — p) 

P = o 


The function x ) which will generate an orthonormal basis is obtained via 

2m — 1 

^i x ) = Y 9p+M 2x - V ) 

p = o 
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with g p = (— l) p 1 h 2 m- P +i, P = 1, . . . ,2m and / <p(x)dx = 1. 
The coefficients {/i p }p=j are generally chosen so that 


ipj ik (x) = 2 *ip ( 2 J x-A:+1), 

for j,k integers, form an orthonormal basis and in addition, the function x ) has m van- 

ishing moments 

J Tp(x)x e dx = 0, i = 0, 1, . . . ,m — 1. 

Also we define 

ifjk = 2"s<^(2 - - , x — k + 1). 

Finally, we assume that there exists a real constant r m (ri = |) such that the following 
conditions are satisfied: 


J <p(x + T m )x e dx = 0 for £=l,...,m— 1, 


and / (p(x)dx = 1. 

In this case the quadrature formula becomes: 


1 


S k ~ /Tt(/( 




k-l+r„ 

N 


•) + 0(iV (m+1) )) 


and the initial discretization error is 0(A r- b 7l+1 )) up to uniform translation. 

The decomposition of the vector {sj, . . . , s^} into the basis we use to compute with 
comes via 

K}— > {4 } — «} ♦ K) 

\ K}\ {4} -\ K}- 

This is implemented in O(iV) operations using: 

p—2rn 

4 = Z * P 4+2fc-i 

p=l 

p—2m 

4 = Z 5p4+2Jt-i 

p=i 

and the s k , d k are viewed as periodic sequences with period 2 V ~K 
The orthonormal basis consists of 

2 4 

The inverse mapping can also be done in 0(N ) operations. 
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Each of the s 3 k is thought of as approximating 

4 = J f(x)<pjk{x)dx = 

2-(^>[/(2- l ' +i (fc-l+r m )) 

_!_£)( jy(-"+iK m + 1 ) )j 

while each d J k is thought of as approximating 

4 = J f{x)^ jk (x)dx. 

The numerical procedure effectively transforms the approximate discretization of the 
matrix G(xj,yk,t n ) which is ( A n )j k . Estimate (2.8) (corresponding to (4.5) and (4.6) of [1], 
uniform in all parameters, indicates (via an argument of [1]) that truncating A n by removing 
elements of a band of width b > 2 m around a shifted diagonal (and its periodic extension) 
i.e., those for which 

| j — k — aAn| > b > 2m, 
which replaces A n by A n ' b , leads to an estimate 

|M n -A*‘||<£tog(J'0 

for C depending only on G. 

It also follows easily that for large N and fixed precision e, only 0(N log N) elements 
will be greater than e. Alternatively, by discarding all elements that are smaller than a fixed 
threshhold we compress it to 0(N\og N) elements. Again following the discussion in [1], we 
note that this naive approach is to construct the full matrix in the wavelet basis and then 
to threshhold. Clearly this is an 0(N 2 ) operation. 

Since we have, a priori, the structure of the singularities of the matrix A v the relevant 
coefficients can be evaluated by using the quadrature formulas. Estimate (2.8) guarantees 
that this procedure requires 0(N log N) operations. 


Remark R2. It is interesting to note that so called unstable difference schemes can be used 
without any drastic loss of efficiency. If (2.1) is approximated by, 

uf 1 = uj - - uJ_,)/2, 


(2.9) 


Uj — Uo(xj), j — 1,2, . . . , N 


the algorithm is not stable for any fixed A > 0, see e.g. [7]. 
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The approximation does converge if At < CAx 2 , (A < CAx ) with an amplification 
factor 1 + O(At). The number of timesteps for t = 0( 1) calculation will be large, n = 
0(Ax~ 2 ) = 0(N 2 ). This is devastating for the standard explicit algorithm (1.2) but will 
only affect the complexity of (1.4) by a constant factor. The number of iterations (m in 
(1.4)) will increase from log(iV) to log(iV 2 ). 

Our approach is in general not as favorable for multidimensional hyperbolic systems, 

d 

d t u + Aj(x)d X} u = f(x), x € R d , 

(2.10) i =1 

u(x , 0) = uo(x). 

When u is a scalar or if the system can be diagonalized the algorithm (1.4) works well. The 
solution is given by integration along characteristics and the support of the Green’s function 
is a small number of points (see Remark (Rl) above). In the idealized case each row of 
A v consists of a fixed number of 6-functions. Its wavelet representation will have log(A^ d ) 
nonzero terms. The overall complexity for (1.4) is then 0((log N) 3 N d ) when the knowledge 
of the location of the 6-functions is used. This is better than the standard 0(N d+1 ) estimate. 

In general, however, the Green’s function for (2.6) has a support with positive volume 
in R d and with a singular support of positive measure in Hausdorff dimension d — 1. The 
representation of the singular support consists of 0(A r<1-1 )6-functions in each row of A". 
This corresponds to C?(log(A r )A r<i_1 ) wavelets and the overall algorithm contains at least 
(0(logiV) 2 ./V 2 ‘ i-1 ) wavelets. 

For general multidimensional problems the new algorithm is still of interest in special 
cases, e.g., if the solution is needed only at a fixed number of points and if it is needed for a 
large number of different data uq, /. 


3. Parabolic Problems 

The Green’s function for parabolic problems is smooth in contrast to the hyperbolic case. 
The pure initial value problem for the heat equation, 

d t u = Au, t > 0, x € R d , 


(3.1) 


u(x,0) = U 0 (x), 


has a solution of the form, 

(3.2) u(x,t) = (M)~ d/2 [ exp(-\x - y\ 2 / it)u 0 (y)dy . 

JRd 
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In bounded domains the kernel has to be changed slightly depending on the boundary 
conditions. For positive t(— nAt) each row in A n is always an approximation of segments of 
regular functions. 

Our new technique is in general more favorable for parabolic problems than hyperbolic 
ones. The structure of the matrix B in (1.4) is simpler. When t increases the kernel becomes 
smoother and ajk can be truncated to zero for all k when j is large enough. 

Explicit methods for (3.1) also requires more operations than for hyperbolic problems 
when the standard method is used. This follows from the parabolic stability requirement, 

(3.3) At < const. |Ax| 2 . 

The new technique is only marginally affected by the constraint (3.3). Compare here the 
discussion above for unstable hyperbolic methods. 

In more general higher order multidimensional parabolic cases the fundamental solution 
of, e.g., 

u t + (— A) d u = 0 
is 

1 t CO 

G d (x, t) = — f d£ exp (i£ • x - \(\ 2 d t). 

Zll J — OO 

This is merely a multidimensional and rescaled version of the fundamental solution used in 
(2.8), and a simpler, but multidimensional version of (2.8) is just: 

\\x\ m+ 1 D?G d (x,t)\<C md . 

Moreover C md is arbitrarily small if t is large enough (this of course requires the nonexistence 
or other special behavior of lower order terms). 

The matrix compression technique is easy here (for periodic problems without boundary 
conditions) because the significant terms of [A v \ lie near the main diagonal and its periodic 
extension in one dimension. In two space dimensions (as is usual for elliptic operators), we 
also need to consider diagonals i = j ± kN for 0 < k < d. Recall A is an N 2 x N 2 matrix in 
2D. 

It is clear that a priori thresholding (to obtain 0(e) precision) near the image of these 
diagonals will give us an 0(N d (\og N) 3 ) operation for each evaluation of the solution, where 
d is the number of space dimensions for the problem. 

4. Numerical Experiments 

The algorithm (1.4) was applied to hyperbolic problems in one space dimensions and to 
one and two dimensional parabolic problems. Various difference approximations and wavelet 
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spaces were used. We present results concerning the accuracy of the calculations and the 
sparsity of ( SAS -1 ) n . 


4.1 Hyperbolic problems. Consider the following scalar hyperbolic problem: 

d t u + a{x)d x u = /(x) 


(4.1a) 


u(x,0) = u 0 (x) 

with periodic boundary conditions (0 < x < 1). We made the following choices: 
(4.1b) a(x) = 0.5 + 0.115 sin(47rx) 


(4.1c) 


f(x) = COs(47Tx) 


(4.1d) u o(x) = sin(4xx). 

In the discretization, Ax = 1/1024 and At /Ax = 1. The wavelet transform operator 
S uses the Daubechies-8 wavelets, which have 8 coefficients and have 4 vanishing moments. 
Finite difference schemes of order 1,2, 3, 4, and 5 of accuracy are tested. 

These finite difference schemes are obtained as follows. In each interval 

(4.2) i - {x/(i / — l)Ax < x < i/Ax} 

a polynomial of degree k is constructed. This polynomial interpolates the two points 
(x„_i ,<_!) and {x v ,ul) and k - 1 of its neighbors. If k is even these interpolation points 
go from x„_k to x ll+ k. If k is odd they go from to x^*^. This gives us a 

reconstruction function which is a polynomial of degree k in each /„ i and is continuous, 
but generally not differentiable at the boundary points x u -\ and x„. We call this function 
R n ' k (x) 

To approximate (4.1) at the grid points (x 1/ ,/ n+1 ) we solve (4.1) “exactly” with initial 
data 

(4.3) u Ax {x,t n ) = R n ' k {x) 

for t n < t < t n+1 , evaluate the solution at (x v ,f n+1 ), and set u" +1 = ua^x^,^ 1 ). We 
require ; go the solution depends only on data in i if a(x) > 0 and if 

a(x) < 0. 
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In the special case when a(x) = a, constant, then 


(4.4) 


u n+l = Rn,k( Xv _ a A<) 



— s))ds 


In the case when / = 0 we get some familiar schemes: For k = 1 this is just the first order 
accurate upwind difference scheme (2.4). For k — 2 this is just the classical Lax-Wendroff 
second order accurate three point scheme, see e.g. [7]. For k = 3,4,5 the schemes are less 
studied, but are known to be L 2 stable, see e.g. [9] and the references therein. 

For variable coefficients the result is 


(4.5a) 


where x v {t) solves 
(4.5b) 


u Ax (x^t n+1 ) = R n ' k (x v {t n )) 

+ / f(Xv{t n+1 - s))ds 
Jt n 

= a(x„), t n < t < t n+1 


(4.5c) x„(f n+1 ) = x u . 

A fourth order Runge-Kutta method is used to integrate the O.D.E. (4.5b, c) and Simpson’s 
rule is used to evaluate the integral in (4.5a). The result of this approximation to the right 
side of (4.5a) is defined to be u" +1 . 

Returning to the present case the computations ran 13 steps until t = 4, that is, 
(S'/IS' -1 ) 2 ’ 3 was computed. 

At each step n the number of elements of A n and (S'AS' -1 )" whose absolute values are 
greater than 10 -4 is shown in table 1. This is for methods whose order of accuracies go from 
one through five. The results are also plotted on Figure 1. 

These significant elements are located near the sub-diagonal corresponding to the char- 
acteristic curve which is known a priori. The image of these locations in (S AS -1 ) 71 , shown 
on figure 2, has total length of 0(N log N) elements where N = 1024. 

In the computation of (SAS* 1 ) 11 , first, from the knowledge of the PDE, we figure out the 
structure of the singularities of A and its image in (S AS~ l ) n . Then we compute (5'A5' _1 ) 2n = 
(5’j4S ,_1 ) n * (5A5 1-1 )" considering only the elements in a neighborhood of the singularities. 
In particular, we define the neighborhood of a singularity to be locations whose distance 
from the singularity are less than or equal to 5. If the singularities lie on a subdiagonal and 
its periodic extension its neighborhood form a subband of bandwidth 11 (the wavelet filters 
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have 8 elements). This bandwidth is independent of the time t (the step n) and the size of 
the problem. The errors due to the subband truncation, measured by ||u n — u n ||/|[u n ||, are 
shown in table 2b. Table 2a shows the relative error between the subband truncation and 
the exact solution. Here and throughout, “|| • ||” denotes the £ 2 norm. Table 2c shows the 
relative error between the subband truncation and untruncated under grid refinement for the 
various orders. Unsurprisingly, since the relative length of the subband which is preserved 
decreases linearly with grid size, the error increases, but only slightly under this process. 

We note that the compression (as seen in Figure 1 and Table 1) is better for odd order 
than for even order schemes. This is perhaps not surprising since (2.7) models schemes of 
odd order accuracy. Singularities behave a bit differently for even order (say order = 2p) 
schemes. These are modeled by 


(4.6) 


u t + au x 


= £ P { Ax) 2p 



u 


+(-iyk p (Ax) 2 > +1 



u 


where k p > 0 and l p are nonzero constants. The odd order dispersive term above may tend 
to spread singularities of the fundamental solution spuriously. 

Finally table 3 shows the relative error due to truncation when the band width of the 
subband is 9, 11, and 13 for the methods of first and second order. Figures 3a and 3b 
compare the truncated versus the approximate solutions due to truncation of bandwidth 9 
for the first and second order methods (the truncated graphs are dotted). 


4.2 Unstable Schemes. For theoretical interest, we apply the method to a finite difference 
scheme which is unstable for -^ = A > 0 


(4.7a) 

= a” - A(u” +1 - u”_i)/2, 

(4.7b) 

u] = U 0 (Xj). 


The amplification factor of this scheme is 
(4.8) 1 — \i sin0 = r(e ,e ), — 7r < 0 < 1 


so 

This means that if 
(4.9) 

13 


|r(e‘ fl )| = (1 + A 2 sin 2 0)K 
A< < 2c(Ax) 2 


for some c > 0, then 
(4.10) 


||/l n ||,i < c" ,Al . 


The restriction (4.9) means that the operation count for this explicit method would be 
0(N 3 ) if we were silly enough to use it. However our compression method allows for an 
operation count of 0(N (log N) 3 ) for the reasons described above. 

Table 4 shows the number of elements in A n and (SAS -1 )" whose absolute values are 
greater than 10 -3 . We choose a bigger threshold here since we took = 1 and nAt = 2, 
so ||A"||, as estimated in (4.10) grows to be roughly 10 when we are finished computing. 

The error as measured by N (subband truncation using bandwidth 11) was 0.0136. 

We also performed convergence studies as we refined the grid for this method. Figures 
(4a,b,c) compare the numerical (untruncated) using dots versus exact solution for m = 
128, 256,512 grid points. The result indicates a second order method, as it should, since At — 
(Ax) 2 . Figures (5a,b,c) compare the truncated bandwidth (using dots) vs the untruncated 
for this method for m ^ 128,256, and 512 grid points. 

The relative error decreases with mesh refinement. The truncation error equation associ- 
ated with this scheme involves limited antidiffusion. Perhaps this accounts for this behavior. 


4.3 System of Hyperbolic Equations. We apply the method to solving the system of 
hyperbolic equations: 


(4.11a) 

on 0 < x < 1, t > 0 with 


(4.11b) 



+ 


a 0 
0 —a 


d x 


v 

w 


0 

0 


the boundary conditions and initial conditions: 


w(0,f) = tn(0,f) 


MM) = KM) 

v(x,0) = v 0 ( x ) 


u>(x,0) — u>o(x) 

the coefficient a is chosen to be constant: 


a = 0.115. 


The numerical method used is the first order accurate upwind method described above. 
The results are similar to the scalar case, except the structure of the singularities in the 
matrices is more complicated. We have to keep track of reflections of singularities at the 
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boundaries which is quite simple in this case. The number of elements in A n and (5v45' _1 )’ 1 
whose absolute values are greater than 10 -4 is shown on table 5, and is plotted on figure 6. 
The relative error due to the subband of width 11 truncation, measured by ||u n — u n ||/||u n ||, 

is 0.0149. 

The structure of the elements whose absolute values are greater than 10 -4 of A 2048 and 
(5'i45 ,_1 ) 2048 is shown in figures (7a, c), while Figure (7b) shows the image of a subband of 
bandwidth 11 in (i>/45 -1 ) 2048 . 


4.4 Parabolic Problems. We do experiments on the following parabolic problem: 

d t u — d x (a(x)d x u ) + /(x) 


(4.12) 


u(x,0) = uo(x) 

with periodic boundary conditions (0 < x < 1). We made the following choices: 
a(x) = 0.5 + 0.25 sin(27rx) 

/(x) = —7 r 2 cos(27rx) 2 + 7 t 2 (0.5 + 0.25sin(27rx)) sin(27rx) 


uq(x) = sin(47rx). 


The discrete setting and the wavelets are the same as in the hyperbolic problem. We use 
the simple explicit central difference scheme (4.13) 

u ” +1 = u n j + . . . 2 A_ (a(xj) A + u j) 

(4.13) 

+ At/(x_,) 


where 

A^uj = T(« iTl - uj ) 

with Af/(Ax) 2 = 0.25. The number of significant elements in A n and (Sb45' -1 ) n is shown on 
table 6, and is plotted on figure 8. 

For the parabolic problem, the large elements of A are in the neighborhood of the main 
diagonal. Their wavelet transform image is shown in figure 9. The relative error due to 
subband truncation was 0.0025. 


4.5 Two-dimensional Parabolic Problems. We consider the following problem: 

diZl — T “f 


u(x,y,0) = u 0 (x,y) 
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with periodic boundary conditions (0 < x < 1, 0<y<l). We choose 

a u (x, y) — 0.5 + 0.25 sin(27rx) 

012 (x,y) = 0.115sin(27rx) cos(27ry) 

— 0.5 + 0.25 cos(27ry) 
uo(x,y ) = sin(47rx) + cos(87rx). 

We use a standard two-dimensional explicit central difference scheme. The two-dimensional 
data Uj'k, j = 1 ... Ni, k = l ... N? forms a one- dimensional vector in the following way 

. . . ui tN2 , tt 2 , i . . . u 2 ,n 2 , • • ■ , Un u i . . . u Nl )N2 } . 

To reduce the size of the problem, N 2 is much less than jVj. In particular we took Ni = 
128, N 2 = 8 that is, Ax = Ay = |. 

The compression worked quite well. Table 7 shows the number of elements in A n on 
(SAS -1 )" whose absolute values are greater than 10 -4 . The relative error due to subband 
truncation was 0.0066. 
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Table 1: Hyperbolic equation: the number of elements in A * and (SAS~*) n 
whose absolute values are greater than 10" 4 
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Table 2: Hyperbolic equation the err ora , meaiured by ||u* - ti*||/||ti*||, (a) compare with 
the eact aolution; (b) due to the truncation only: (c) due to the truncation only under grid 
refinement. 
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Table 3; Errors measured by jl due to truncation for various band widths and first 

and second order. 
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Table 4: Hyperbolic equation ‘"unstable scheme*; the number of elements in A " and 
(5-45* 3 ) n whose absolute values are greater than 10"* 


.HI 

ma 


§■] 

T'TTS 


HU 

■Ml 

■iga 

^■3 

KkEl 


HD 

M3EM 

Kina 

■U 

■Ml 


mm 

caa 


mm 

caa 

WMZMM 

■£1 


■Mga 

■21 

bebi 


ca 

Eaa 


IMl 

BMI 

khlou 


EBEa 

■E£BI 


Table 5: System of bj-perbolic equations: The number of elements in A " and (5A5“ , ) n 
whose absolute values are greater than 10” 4 . 
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Table 6 Parabolic equation: the number of elements in A* and (5A5’ 1 )* whose absolute 
values are greater than 10“ 4 
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Table 7: 2D*parabolic equation the number of elements in A " and ($.4S“ 3 )" whose 
absolute values axe greater than 10“ 4 
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order 1, width 9 



Figure 3a Thincation versus nontruncated approximate solution, first order method trun- 
cated at bandwidth 9 (Truncated is dottedj. 


order 2, width 9 



Figure 3b: Truncated versus nontruncated approximate aolution, aecoDd order method, 
truncated at bandwidth 9. (Truncated is dotted). 
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Figure 6: System of hyperbolic equations: The number at dementi is A m end (SAS~ J ) 9 
whose absolute values are greater than 10~ 4 . 



Figure 7a System of hyperbolic equations patters of significant elemenu'JjalO" 4 lor 
X", n - 2048 ** 



Figure 7b: System of hyperbolic equations pattern of significant element for (SAS~ l )',«■ 
2048. image of bandwidth 11 around singular support 
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Figure 0: Parabolic aquation: tbc pattern of ngnificant tknjeat* in (SAS* 1 )". 
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